{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Математическая статистика\n",
    "## Практикум 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as sps\n",
    "import scipy.stats as stats\n",
    "import ipywidgets as widgets\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "%matplotlib inline\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Задача 1**  \n",
    "Пусть $x_1,\\ldots,x_n$ — реализация из модели сдвига показательного закона\n",
    "\tс плотностью распределения \n",
    "\t$$\n",
    "\t\tf_\\theta(u) = e^{-\\lambda(u-\\theta)}\\cdot I_{\\{u\\geq\\theta\\}} \n",
    "\t\t= \n",
    "\t\t\\begin{cases}\n",
    "\t\t \te^{-\\lambda(u-\\theta)}, & u\\geq\\theta,\\\\\n",
    "\t\t\t0, & u<\\theta.\n",
    "\t\t\\end{cases}\n",
    "\t$$\t\n",
    "Оценить $\\theta$ с помощью метода моментов и метода максимального правдоподобия. Реализуйте эту задачу в Python:\n",
    "\n",
    "* сгенерируйте $\\theta$ из равномерного распределения на  $[1,N+1]$ и возьмите $\\lambda=|N-10|/N$, где `N` - ваш  номер в журнале группы и ФИО(нумерация начинается с 01,02,03... и т.д.)\n",
    "* сгенерируйте выборку из распределения с плотностью $f_\\theta(u)$ размера  $n=1\\,0000$; \n",
    "* постройте алгоритм нахождения оценки методом моментов и методом правдоподобия (мы такой пример решали на занятии), а потом найдите значения оценок метода моментов и метода максимального правдоподобия, которые генерируются автоматически Питоном."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(100) # зафиксируем seed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Seed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для того, чтобы каждый раз при использовании генератора псевдослучайных чисел получать идентичные последовательности, используется функция set.seed (от set - задать, установить, и seed - начальное число). Как следует из названия, эта функция фиксирует число, служащее начальной точкой для запуска алгоритма генерации псевдослучайных чисел. В качестве аргумента функции указывают любое целое число (не важно, какое именно). \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **np.random.seed(<число>)** —\n",
    "настраивает генератор случайных чисел на новую последовательность. Если данная функция в программе не вызывается, по умолчанию используется системное время. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В качестве примера приведем еще фукции для распределений, которые мы упоминали в курсе по теории вероятностей:\n",
    "* **np.random.binomial(n, p[, size])** — биномиальное распределение  \n",
    "* **np.random.exponential([scale, size])** — экспоненциальное распределение  \n",
    "* **np.random.poisson([lam, size])** — распределение Пуассона  \n",
    "* **np.random.standard_cauchy([size])** — распределение Коши  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Обратите внимание на то, что в np.random.exponential параметр scale = 1/$\\lambda$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(14, 0.2857142857142857)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 14 # номер варианта\n",
    "𝜆 = np.abs(N - 10) / N\n",
    "\n",
    "N, 𝜆"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10000 # размер выборки\n",
    "scale = 1 / 𝜆\n",
    "θ = np.random.uniform(1, N + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = θ + np.random.exponential(scale, size = n) # генерируем выборку размера n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "esimate_mme = np.mean(sample - scale) # оценка методом моментов\n",
    "esimate_mle = np.min(sample) # оценка методом максимального правдоподобия"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Истинное значение параметра: 8.607669185073515\n",
      "Оценка метода моментов: 8.59730674766141\n",
      "Оценка метода максимального правдоподобия: 8.607720350439438\n"
     ]
    }
   ],
   "source": [
    "print(f'Истинное значение параметра: {θ}')\n",
    "print(f'Оценка метода моментов: {esimate_mme}')\n",
    "print(f'Оценка метода максимального правдоподобия: {esimate_mle}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Задача 2**  \n",
    "Мы нашли точечные оценки методом моментов и методом правдоподобия, теперь надо оценить используя соответствующий параметрическому случаю критерий согласия, согласуется ли полученные распределения с данной выборкой. Good-to-fit критерии или критерии согласия с простой гипотезой. Используйте любой критерий, который проверит является ли распределение экспоненциальным с заданным показателем. Научитесь подбирать подходящие критерии. Какие критерии бывают можно посмотреть в книжке Медведева, Ивченко."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Результат с оценкой метода моментов\n",
      "Статистика: 0.0079\n",
      "p-value: 0.9139260465470475\n",
      "\n",
      "Результат с оценкой метода максимального правдоподобия\n",
      "Статистика: 0.0088\n",
      "p-value = 0.8335435086261689\n"
     ]
    }
   ],
   "source": [
    "# оценка метода моментов\n",
    "mme = esimate_mme + np.random.exponential(scale, size = n)\n",
    "kstest_mme = sps.ks_2samp(sample, mme)\n",
    "print('Результат с оценкой метода моментов')\n",
    "print(f'Статистика: {kstest_mme.statistic}')\n",
    "print(f'p-value: {kstest_mme.pvalue}\\n')\n",
    "\n",
    "# оценка метода максимального правдоподобия\n",
    "mle = esimate_mle + np.random.exponential(scale, size = n)\n",
    "kstest_mle = sps.ks_2samp(sample, mle)\n",
    "print('Результат с оценкой метода максимального правдоподобия')\n",
    "print(f'Статистика: {kstest_mle.statistic}')\n",
    "print(f'p-value = {kstest_mle.pvalue}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$p-value > 0.05 ⇒$ нет оснований для отклонении гипотезы о том, что выборки принадлежат одному типу распределения.\n",
    "\n",
    "Оценки методами 'моментов' и 'максимального правдоподобия' дают согласованные результаты."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Задача 3** \n",
    "Сравнить на выборках размера 100 для $\\mathcal{N}(\\theta,4)$ доверительные интервалы:\n",
    "(1) теоретический, (2) на основе параметрического бутстрэпа, (3) на основе непараметрического бутстрэпа. Сам параметр $\\theta$ сгенерировать из равномерного распределения на $[-5,5]$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100 # размер выборки \n",
    "α = 0.05 # уровень значимости\n",
    "\n",
    "θ = np.random.uniform(-5, 5)\n",
    "σ = np.sqrt(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = np.random.normal(θ, σ, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Значение θ: -3.6100308551787066\n"
     ]
    }
   ],
   "source": [
    "print(f'Значение θ: {θ}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Теоретический доверительный интервал"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Напомним, что теоретический доверительный интервал вычисляется следующим образом: \n",
    "\n",
    "$$\n",
    "\\mathbb{P}\\left( \\bar{X} - \\frac{c_{1-\\alpha/2}\\sigma}{\\sqrt{n}} < \\mu < \\bar{X} + \\frac{c_{1-\\alpha/2}\\sigma}{\\sqrt{n}} \\right) = 1-\\alpha,\n",
    "$$\n",
    "где $c_{\\alpha}$ — квантиль распределения $\\mathcal{N}(0,1)$ уровня $\\alpha$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Вычисляем теоретический доверительный интервал используя функции $stats.norm.ppf(1-alpha/2)$ для вычисления кваниля $c_{\\alpha}$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Теоретический доверительный интервал: (-3.888687526547435, -3.0792696230134338)\n"
     ]
    }
   ],
   "source": [
    "n = len(sample)\n",
    "avg = np.mean(sample)\n",
    "std = np.std(sample, ddof = 1)\n",
    "\n",
    "quantile = stats.norm.ppf(1 - α / 2)\n",
    "\n",
    "confidence_interval = (avg - quantile * std / np.sqrt(n), avg + quantile * std / np.sqrt(n))\n",
    "print(f'Теоретический доверительный интервал: {confidence_interval}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Доверительный интервал на основе параметрического бутстрэпа"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_of_bootstrap_samples = 10 # количество бутстрэп-выборок\n",
    "size_of_bootstrap_samples = 20 # размер бутстрэп-выборок"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Параметрический бутстрэп доверительный интервал: [-4.140715248471226, -2.8949365807765646]\n"
     ]
    }
   ],
   "source": [
    "bootstrap_samples = np.random.normal(np.mean(sample), σ, size = [number_of_bootstrap_samples, size_of_bootstrap_samples]) \n",
    "bootstrap_estimates = np.apply_along_axis(np.mean, 1, bootstrap_samples)\n",
    "\n",
    "parametric = [np.quantile(bootstrap_estimates, α / 2), np.quantile(bootstrap_estimates, 1 - α / 2)]\n",
    "\n",
    "print(f'Параметрический бутстрэп доверительный интервал: {parametric}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Доверительный интервал на основе непараметрического бутстрэпа"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_of_bootstrap_samples = 10 # количество бутстрэп-выборок\n",
    "size_of_bootstrap_samples = 20 # размер бутстрэп-выборок"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Непараметрический бутстрэп доверительный интервал: [-4.634893425709553, -3.0533332066753043]\n"
     ]
    }
   ],
   "source": [
    "bootstrap_samples = np.random.choice(sample, size = [number_of_bootstrap_samples, size_of_bootstrap_samples]) \n",
    "bootstrap_estimates = np.apply_along_axis(np.mean, 1, bootstrap_samples)\n",
    "\n",
    "nonparametric = [np.quantile(bootstrap_estimates, α / 2), np.quantile(bootstrap_estimates, 1 - α / 2)]\n",
    "print(\"Непараметрический бутстрэп доверительный интервал:\", nonparametric)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Выводы: как сравнить полученные доверительные интервалы? Как сравнить методы? Какие лучше и почему?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Значение θ: -3.6100308551787066\n",
      "\n",
      "Теоретический доверительный интервал: (-3.888687526547435, -3.0792696230134338)\n",
      "Параметрический бутстрэп доверительный интервал: [-4.140715248471226, -2.8949365807765646]\n",
      "Непараметрический бутстрэп доверительный интервал: [-4.634893425709553, -3.0533332066753043]\n"
     ]
    }
   ],
   "source": [
    "print(f'Значение θ: {θ}\\n')\n",
    "print(f'Теоретический доверительный интервал: {confidence_interval}')\n",
    "print(f'Параметрический бутстрэп доверительный интервал: {parametric}')\n",
    "print(f'Непараметрический бутстрэп доверительный интервал: {nonparametric}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<strong>Как сравнить полученные доверительные интервалы?</strong>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сравниваем интервалы по ширине"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<strong>Как сравнить методы?</strong>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Доврительные интервалы содержат значение параметра θ. Разница в ширине:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Значение θ: -3.6100308551787066\n",
      "\n",
      "Теоретический доверительный интервал: (-3.888687526547435, -3.0792696230134338), θ входит\n",
      " - Ширина 0.8094179035340012\n",
      "\n",
      "Параметрический бутстрэп доверительный интервал: [-4.140715248471226, -2.8949365807765646], θ входит\n",
      " - Ширина 1.2457786676946614\n",
      "\n",
      "Непараметрический бутстрэп доверительный интервал: [-4.634893425709553, -3.0533332066753043], θ входит\n",
      " - Ширина 1.5815602190342486\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(f'Значение θ: {θ}\\n')\n",
    "\n",
    "def is_inside(interval):\n",
    "    return 'θ ' + ('' if interval[0] <= θ <= interval[1] else 'не ') + 'входит'\n",
    "\n",
    "def print_stat(message, interval):\n",
    "    print(f'{message}: {interval}, {is_inside(interval)}')\n",
    "    print(f' - Ширина {interval[1] - interval[0]}\\n')\n",
    "\n",
    "print_stat('Теоретический доверительный интервал', confidence_interval)\n",
    "print_stat('Параметрический бутстрэп доверительный интервал', parametric)\n",
    "print_stat('Непараметрический бутстрэп доверительный интервал', nonparametric)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<strong>Какие лучше и почему?</strong>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "θ входит во все интервалы, а теоретический самый узкий - поэтому он подходит лучше всего."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
